#####################################
# Application
# Government formation
# Q1
#####################################
rm(list=ls())
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
source("cbq_function.R")
    
# load data
dat0 <- read.table("formation_new.tab",sep = "	",header=T)
dat <- dat0[order(dat0$case,dat0$realg),]
y <- dat$realg
Y <- ifelse(y==0,-1,1)
N <- length(y)
ind <- as.integer(as.factor(dat$case))
n_grp <- max(ind)
X <- as.matrix(subset(dat,select = c(minor,
    minwin,
    numpar,
    dompar,
    median,
    gdiv1,
    mgodiv1,
    prevpm,
    sq,
    mginvest,
    anmax2,
    pro,
    anti)))
n_coef <- dim(X)[2]
dat$ind <- ind
    
####################################
# estimate with CBQ
ms_cbq_stan <- cbq(N = N,
    n_covariate = n_coef,
    X = X,
    Y = Y,
    N_indx = n_grp,
    ind = ind,
    qtl = 1,
    nchain = 5,
    niter = 1000,
    thin = 1)
    
save(ms_cbq_stan,file="cbq_q1.RData")
